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Abstract 

Mutations are central to evolution, providing the genetic variation upon which selection acts. A mutation's effect on the 
suitability of a gene to perform a particular function (gene fitness) can be positive, negative, or neutral. Knowledge of the 
distribution of fitness effects (DFE) of mutations is fundamental for understanding evolutionary dynamics, molecular* 
level genetic variation, complex genetic disease, the accumulation of deleterious mutations, and the molecular clock. We 
present comprehensive DFEs for point and codon mutants of the Escherichia coli TEM-1 |$-lactamase gene and missense 
mutations in the TEM-1 protein. These DFEs provide insight into the inherent benefits of the genetic code's architecture, 
support for the hypothesis that mRNA stability dictates codon usage at the beginning of genes, an extensive framework 
for understanding protein mutational tolerance, and evidence that mutational effects on protein thermodynamic sta- 
bility shape the DFE. Contrary to prevailing expectations, we find that deleterious effects of mutation primarily arise from 
a decrease in specific protein activity and not cellular protein levels. 
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Introduction 

The fitness landscape model for protein evolution, as first 
conceptualized by Smith in 1970 (Smith 1970) and general- 
ized by others (Orr 2005), imagines evolution as a process by 
which a sequence moves by stochastic processes from its 
wild-type sequence through fitter and fitter sequences until 
the sequence reaches a local fitness optimum. The nature of 
the fitness landscape determines the dynamics of evolution 
and fundamentally shapes what is and is not possible in evo- 
lution. Much has been learned from theoretical studies and 
small-scale interrogations of real fitness landscapes (Eyre- 
Walker and Keightley 2007). However, we still lack a system- 
atic, assumption-free, experimental determination of the dis- 
tribution of fitness effects (DFE) for all mutations of a gene 
performing its native function in its native host. The situation 
is akin to having a small set of aerial photographs of a geo- 
graphical area versus having comprehensive satellite coverage 
such as provided by Google Earth. 

We sought to provide a comprehensive, quantitative 
description of a fitness landscape corresponding to a gene 
and its nearest neighbors in both DNA and protein sequence 
space (i.e., the set of all sequences that differ by a single bp, 
codon, or amino acid) but avoid or ameliorate current limi- 
tations of large-scale measurements of fitness. Growth com- 
petition experiments or experiments in which alleles are 
enriched based on a threshold for function are the current 
state of the art (Fowler et al. 201 0; Araya et al. 201 2; Deng et al. 
2012; McLaughlin et al. 2012; Schlinkmann et al. 2012; 
Whitehead et al. 2012; Roscoe et al. 2013; Starita et al. 
2013). To varying degrees such experiments offer a direct 
"head-to-head" comparison of alleles but suffer four 



significant limitations. First, most studies utilize nonnative 
reporter assays (e.g., phage display, cell surface display, and 
two-hybrid systems) in which the gene or gene fragment is 
removed from its native context and host and fused to an- 
other gene (but see Deng et al. 2012; Roscoe et al. 2013). 
Second, population size can affect the measured value of fit- 
ness due to stochastic effects. Third, these experiments have 
limited ability to measure fitness for low fitness alleles because 
such alleles are depleted during the course of the experiment. 
For example, Roscoe et al. (2013) were unable to reproducibly 
measure the fitness of ubiquitin point mutants with a fitness 
below approximately 40% of wild-type due to their rapid 
depletion in growth competition experiments. Thus, al- 
though such experiments tell us the location of valleys in 
the landscape, they cannot tell us anything about what the 
valleys look like. Fourth, the fitness measurements are subject 
to the extent and underlying form of genotype-by-environ- 
ment interactions. For example, the fitness of an antibiotic 
resistance gene measured by a growth competition experi- 
ment will be a function of the arbitrary selective pressure used 
in the experiment (the antibiotic concentration). Alleles con- 
ferring resistance far above or below the level necessary for 
growth at one antibiotic concentration may show no fitness 
difference in that environment yet show significant differ- 
ences at a different antibiotic concentration closer to their 
resistance limit. We desired to decouple fitness from geno- 
type-by-environment interactions as much as possible to 
quantify the underlying landscape and thus better under- 
stand a gene's intrinsic evolutionary potential and limitations. 

A key determinant of the fitness landscape is what we are 
defining here as "gene fitness." "Fitness" in the traditional 
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biological definition refers to the extent to which an organism 
is adapted to or able to produce offspring in a particular 
environment (often measured as growth rate in the fitness 
landscape of bacteria). Instead, we are referring to "gene fit- 
ness" as a phenotypic signature — the suitability of a gene to 
provide a particular function. Thus, a gene fitness landscape 
might also be referred to as a phenotypic landscape. Although 
to a large extent this landscape reflects a protein functional 
landscape, we use the term "gene fitness" because our land- 
scape encompasses mutational effects at the DNA, RNA and 
protein level and is a holistic metric of the ability of an allele to 
provide a particular function. However, when we average the 
gene fitness values of synonymous alleles to examine the 
effect of missense mutations, we refer to this as a protein 
fitness landscape. "Protein fitness" is similarly defined as the 
suitability of a protein to provide a particular function. Unless 
otherwise specified, use of the word fitness in relation to a 
specific allele refers to "gene fitness," whereas its use in rela- 
tion to a specific protein refers to "protein fitness." These 
landscapes define the relationship between DNA/protein 
sequence and biological function. The more that function is 
related to growth rate in particular environment, the more 
direct the relationship between gene fitness and organismal 
fitness. For example, in the case of antibiotic resistance genes, 
gene/protein fitness (as measured by minimum inhibitory 
concentration [MIC] of the antibiotic), and organismal fitness 
(as measured by growth rate) correlate to some extent (Deris 
et al. 2013; Jacquier et al. 2013), provided the antibiotic con- 
centration is at the appropriate level for comparing the alleles. 
Here, we present a comprehensive gene fitness landscape for 
the TEA/1-7 (3-lactamase gene and a protein fitness landscape 
for the TEM-1 protein. 

Results and Discussion 

Fitness Landscape of TEM-1 P-Lactamase 
We chose to measure the DFE of the TEM-1 (3-lactamase gene, 
a convenient model for the study of evolution and the fitness 
effects of mutations (Salverda et al. 2010; Soskine and Tawfik 
2010). TEA/1-1 is native to Escherichia coli as a plasmid-borne 
gene (Medeiros 1984), and we examined TEM-1 in this con- 
text. TEM-7 confers high resistance to penicillin antibiotics 
such as ampicillin (Amp). Thus, when E. coli cells bearing 
TEM-1 are challenged to grow in the presence of Amp, alleles 
conferring an enhanced ability to degrade the antibiotic will 
enrich. Thus, Amp resistance is a key determinant of organ- 
ismal fitness in the presence of Amp (Bershtein et al. 2006; 
Weinreich et al. 2006; Jacquier et al. 2013), although assessing 
TEM-1 fitness by measuring Amp resistance does not capture 
organismal fitness differences not associated with antibiotic 
resistance. Previous partial characterizations of the DFE of 
TEM-1 (Bershtein et al. 2006; Deng et al. 2012; Jacquier et al. 
2013) suffer several significant limitations. These studies did 
not characterize the relationship between sequence and fit- 
ness (Bershtein et al. 2006), used error-prone PCR to generate 
mutations that were heavily biased to A/T to C/C transitions 
(80%) (Bershtein et al. 2006; Jacquier et al. 2013), and/or fo- 
cused on the characterization of high fitness alleles with more 



than one mutation and assumed additivity for predicting the 
effect of the individual mutations (Deng et al. 2012). In addi- 
tion, fitness was either measured using either growth compe- 
tition experiments (Deng et al. 2012), which suffer from 
limitations described in the Introduction, or in the coarse- 
grained manner of a MIC assay (Bershtein et al. 2006; Jacquier 
et al. 2013). MIC assays suffer the drawbacks of being low- 
throughput and low-resolution. Alleles with known muta- 
tions must be isolated and tested individually, and MICs are 
measured in discrete values (typically 2-fold increments). For 
example, the resolution of the MIC assay in a study of the 
amoxicillin resistance effects of 18% of the possible amino 
acid substitutions in TEM-1 was insufficient to capture the 
effects of synonymous mutations or to identify any beneficial 
mutations (Jacquier et al. 2013), both of which we readily 
achieve. 

Here, we describe a synthetic biology approach to quantify 
fitness of TEM-7 in a single experiment that avoids or ame- 
liorates the limitations of growth competition experiments 
and MIC assays and allows a comprehensive analysis of the 
DFE. A synthetic biology approach is by definition artificial 
in at least some aspects, but unlike several previous studies 
we measure the DFE of the gene in its native host and do 
not employ gene fusions or artificial reporters of fitness. 
Additionally, because TEM-1 increases the Amp resistance 
of E. coli cells over 1,000-fold, the combination of TEM-7 
and Amp afforded the opportunity to determine the DFE 
over a wide range of fitness values. Our approach decouples 
genotype-by-environment interactions as far as Amp 
resistance is concerned. We quantify TEM-7's underlying fit- 
ness landscape and thus its intrinsic evolutionary potential 
and limitations. This gene fitness landscape is a very signifi- 
cant determinant in the organismal fitness landscape for 
growth of the bacteria in the presence of Amp (Jacquier 
et al. 2013). However, the two types of landscapes are not 
equivalent. 

We determined the DFE for 98.2% (2,536/2,583) of all point 
mutations (i.e., all 1-bp changes) and 83.9% (15,167/18,081) of 
all codon substitutions in the TEM-7 gene (fig. 1). The latter 
includes all 1-, 2-, and 3-bp changes of the 287 codons of 
TEM-1. We also determined the DFE for 95.6% (5,212/5,453) 
of the possible single amino acid substitutions in the corre- 
sponding TEM-1 protein (fig. 2). We excluded insertions and 
deletions (indels) from our analysis of the DFE, an important 
class of mutations that generally have more deleterious effects 
on fitness (Toth-Petroczy and Tawfik 2013). The source of 
TEM-7 variants was a previously described library (CCM-2) 
designed to contain all possible single codon substitutions in 
the TEM-7 gene (i.e., each codon position in the gene could be 
changed to any of the other 63 codons but each allele had 
only one position changed) (Firnberg and Ostermeier 2012). 
To measure gene fitness, we first partitioned the CCM-2 
library into 13 partially overlapping sublibraries based on 
relative Amp resistance using a synthetic gene circuit that 
functions as a tunable bandpass genetic selection for Amp 
resistance (Sohka et al. 2009) (supplementary fig. SI, 
Supplementary Material online). Next, we performed deep 
sequencing on each of the sublibraries, counting how many 
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Fig. 1. Distribution of gene fitness effects (DFE) of mutations in TEM-7. 
(A) The DFE of point mutations (i.e., 1-bp changes in the gene). (B) The 
DFE of all possible codon substitutions (i.e., all 1, 2, and 3 base changes in 
the 287 codons of TEM-7). Gene fitness values for conferring ampicillin 
resistance are presented on a log scale with 0 corresponding to the 
fitness of TEM-1. The contributions of synonymous (red), missense 
(gray), and nonsense (blue) mutations to the DFE are indicated. Gene 
fitness as a function of codon substitution is provided as supplementary 
data SI (Supplementary Material online). 



times each allele appeared in each sublibrary and used these 
statistics to quantity each allele's conferred antibiotic resis- 
tance or gene fitness (w) relative to TEM-1 (supplementary 
fig. S2, Supplementary Material online). We used the fitness 
effects of synonymous mutations to determine an upper limit 
on the error of our fitness measurements (supplementary 
fig. S3 and Materials and Methods, Supplementary Material 
online). Our method enables the accurate fitness quantifica- 
tion of any allele and avoids population size effects because 
the alleles are isolated on a plate. Unlike growth competition 
experiments, the probability of observing an allele in our 
experiment is predominantly independent of its conferred 
fitness. Additionally, the method decouples fitness from 
genotype-by-environment interactions, at least as far as the 
major environmental factor affecting fitness is concerned (i.e., 
the antibiotic concentration). We individually sequenced the 
alleles from 27 randomly selected colonies from two of the 
sublibraries and found these alleles' gene fitness values were in 
the expected range (supplementary fig. S4, Supplementary 
Material online). 



TEM-1's DFE (fig. 1) indicated the gene was fairly robust to 
mutations as nearly one-half (47.3%) of all alleles retained at 
least 50% of the fitness of TEM-1. Among alleles with point 
mutations, 63.8% maintained at least 50% of the fitness of 
TEM-1 (53.2% of the nonsynonymous and 97.2% of the syn- 
onymous point mutations), less than a previous estimation of 
75% (Soskine and Tawfik 2010). Still, a sizable fraction of the 
mutants lost more than 90% of their fitness (19.6% of the 
point mutations and 30.3% of all codon substitutions), 
roughly in line with previous estimates of the frequency of 
mutations having a severe deleterious effect (Camps et al. 
2007). Among point mutants, only 6% of the alleles comple- 
tely lost the ability to provide any Amp resistance and 33% of 
those were nonsense mutations (fig. 1A), which is similar to 
the approximately 8% inactivating mutations found in a pre- 
vious study of an error-prone PCR library of TEM-1 (Soskine 
and Tawfik 2010). Only 7.1% (1,074/15,167) of the alleles and 
7.0% (367/5,212) of the missense mutations increased fitness 
above that of TEM-1 outside the range of the error. The bi- 
modal distribution was qualitatively similar to the DFE of 
randomly chosen point mutations in DNA and RNA viruses 
(Sanjuan et al. 2004; Peris et al. 2010), the DFE of a set of 
induced mutations in yeast (Wloch et al. 2001), the DFE of 
missense mutations of ubiquitin (Roscoe et al. 2013), a sam- 
pling of TEM-1 's DFE for amoxicillin resistance (Jacquier et al. 
2013), and estimations of TEM-1 's DFE for Amp resistance 
(Soskine and Tawfik 2010). 

Benefits of the Genetic Code's Architecture 
TEM-Ts DFE provides evidence that the standard genetic 
code's architecture minimizes the deleterious effects of mu- 
tations and enriches for adaptive mutations. The adaptive 
theory on the origin of the genetic code states that the 
genetic code is arranged to minimize the deleterious effects 
of mutations and mistranslations (Sonneborn 1965; Woese 
1965). This theory predicts that point mutations would be 
less deleterious than 2- or 3-bp substitutions. We have re- 
cently shown this prediction held true for mutations in two 
small genes (HB36 and HB80; Whitehead et al. 2012) that were 
reengineered for a new function in a nonnative organism 
(Fimberg and Ostermeier 2013). Here, we find that this pre- 
diction is also true of a wild-type gene in its native host. The 
median changes in relative gene fitness for 1-, 2-, and 3-bp 
substitutions at a codon position were —0.36, —0.52, and 
—0.63, respectively. More significantly, the frequency of 
point mutations among the alleles with a fitness less than 
0.1 was 35.3% less than that expected if the point muta- 
tions were evenly distributed across all fitness values 
(P= 1.1 x 10~ 40 based on comparison with a hypergeometric 
distribution). For HB36 and HB80, point mutations were 
56.4% and 53.8% depleted from clones with a fitness less 
than 0.1, respectively (P= 1.35 x 10~ 18 and 1.27 x 10~ 21 ) 
(Fimberg and Ostermeier 2013). We interpret this result as 
evidence that the code's arrangement minimizes the fitness 
cost of amino acid substitutions. An alternative explanation is 
that TEM-1, as a product of millions of years of evolution 
under the standard genetic code (Hall and Barlow 2004), 
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Fic. 2. The sequence-function landscape of TEM-1. The heat map indicates the protein fitness values for ampicillin resistance of the indicated amino 
acid substitution. The Ambler consensus numbering system (Ambler et al. 1991) for class A |3 lactamases is used. An asterisk indicates key active site 
residues. For the start codon, fitness values correspond to the average of the codons for the indicated amino acid though methionine is expected to be 
the amino acid incorporated. Protein fitness as a function of missense mutation is provided as supplementary data S2 (Supplementary Material online). 



evolved to minimize the deleterious effects of mutations 
under the rules of this code. 

We also find further support for our hypothesis that the 
standard genetic code's architecture enriches for adaptive 
mutations (Firnberg and Ostermeier 2013). Among the 367 
beneficial missense mutations in TEM-1, 41.1% can be 



achieved by point mutations, 32.5% higher than the 31.0% 
expected if 367 missense mutations were chosen at random 
(P = 8.8x 10~ 6 based on comparison to a hypergeometric 
distribution). Our comprehensive analysis of beneficial muta- 
tions in a natural gene in its native host is the strongest 
evidence yet supporting the hypothesis that the code's 
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arrangement makes adaptive mutations more likely. The role, 
if any, such enrichment played in the origin of the genetic 
code and whether the enrichment is a side effect of the code's 
error minimization bias are difficult questions to answer 
(Fimberg and Ostermeier 2013). 

The Effects of Synonymous Mutations 
The effects of synonymous mutations on protein synthesis 
and fitness have important implications for evolution and 
biotechnology. However, despite an abundance of plausible 
hypotheses, we lack a mechanistic understanding of these 
effects (Plotkin and Kudla 201 1). Our systematic strategy pro- 
vides an assumption-free approach for testing and generating 
these hypotheses. We first examined the gene fitness values of 
725 alleles synonymous to TEM-1. Beneficial and deleterious 
synonymous mutations distributed differently across the 
sequence of TEA/1-7 (fig. 3A). Beneficial mutations occur pri- 
marily in positions 15-30 and 130-260, whereas deleterious 
mutations appeared in clusters in the first half of the gene and 
were almost absent from the second half of the gene. 
No trend in the types of substitutions for either beneficial 
or deleterious effect was apparent other than eight of the ten 
beneficial mutations at Arg codons being to the rare £ coli 
codons AGA (2/10) and AGG (6/10). The pattern of beneficial 
and deleterious synonymous codons indicates the existence 
of regions of TEA/1-7 with suboptimal and less robust mRNA 
properties, respectively. 

We next analyzed the effects of 14,055 synonymous 
substitutions among the set of 15,167 alleles with gene fitness 
measurements (supplementary fig. S3, Supplementary 
Material online). Over the length of the entire gene, CUA 
(Leu), AGG (Arg), and UCG (Ser) provided an average fitness 
advantage over some of their synonymous codons (supple- 
mentary fig. S5, Supplementary Material online), but the ad- 
vantage was only approximately 5%. Interestingly, CUA and 
AGG are rare codons in E. coli. Codon usage often differs in 
the beginning of the gene from the rest of the gene, which has 
been hypothesized to result from a selection against 5' mRNA 
structure and/or a selection for rare codons that provide a 
slower elongation time at the 5' end (Plotkin and Kudla 201 1 ). 
Our data address both these hypotheses. Positions 2-10 in 
TEM-7 had an almost 2-fold broader distribution of synony- 
mous effects compared with any other section of the gene 
(supplementary fig. S3F, Supplementary Material online). 
Within these nine positions, we observed 26-85% mean 
fitness increases for certain codons of Ala, Arg, Gly, Leu, 
Pro, and Ser relative to select synonyms (supplementary fig. 
S6, Supplementary Material online). These synonymous fit- 
ness differences distributed differently among the nine posi- 
tions (supplementary fig. S7, Supplementary Material online). 
Contrary to the slow elongation hypothesis, favored codons 
tended to appear more frequently in the E. coli genome than 
their corresponding disfavored codon (fig. 3B) (Hilterbrand 
et al. 2012). However, none of the 16 observed codon prefer- 
ences were between the most and least frequently used 
codons within a synonym set suggesting that codon usage 
was an inadequate explanation for the observed preferences 
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Fig. 3. Effects of synonymous substitutions. (A) Beneficial and deleteri- 
ous synonymous mutations in TEM-1 are not evenly distributed. Alleles 
synonymous to TEM-1 with a gene fitness significantly higher (red cir- 
cles) or lower (blue squares) than that of TEM-1 are shown. The criteria 
for significance was that the error did not extend into the range fit- 
ness = 1 ±0.1. Error bars provide an upper estimate on error on the 
fitness measurements as described in the text. (B) An analysis of pairs 
of synonymous alleles with mutations in codon positions 2-10 of the 
gene (supplementary figs. S6 and S7, Supplementary Material online) 
revealed that codon preferences tended to be for codons with a higher 
frequency in the Escherichia coli genome. (C) Preferred codons at posi- 
tions 2-10 of the gene were predicted to result in mRNA with less stable 
structures around the initiation codon. Red bar is the mean. **P < 0.01, 
****P < 0.0001 by Student's t-test. 



(e.g., as a result of tRNA abundance). We next calculated the 
folding energy of the mRNA around the initiation codon for 
alleles exhibiting gene fitness differences (Hofacker 2009). In 
almost all cases, favored codons reduced mRNA stability 
around the translation start site compared to disfavored 
codons (fig. 3C). Our findings reinforce recent studies that 
undermine the slow elongation hypothesis (Supek and Smuc 
2010; Charneski and Hurst 2013) and support the theory that 
mRNA structure at the beginning of genes determines the 
translation rate (Bentele et al. 2013; Goodman et al. 2013). 
Like the most recent of these studies (Goodman et al. 2013), 
our study shows how systematic analyses of large synthetic 
libraries is a powerful approach for testing competing 
hypotheses. 

Exceptions to the Standard Genetic Code 
Among the three stop codons, UAG (amber) exhibited non- 
sense suppression (supplementary fig. S8A, Supplementary 
Material online). A 3' flanking purine after the UAG enhanced 
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this suppression (supplementary fig. S8B, Supplementary 
Material online), as has been observed with the amber sup- 
pressor tRNA allele supE (Bossi 1983; Miller and Albertini 
1983). We sequenced the seven tRNAs known to serve as 
amber nonsense suppressors and found that £ coli strain 
SNO301 harbors the supE44 allele, which consists of a dupli- 
cate copy of the glnV tRNA gene, glnX, with the expected 
anticodon mutation (thereby inserting glutamine at UAG 
codons) (Singaravelan et al. 2010). This allele suppressed 
a UAG with a 3' flanking purine at a mean efficiency of 
7-10% (supplementary fig. S8C, Supplementary Material 
online). 

Substitutions for the AUG start codon that provided sig- 
nificant antibiotic resistance (>5% of that of TEM-1) included 
seven of the nine point mutants of AUG (supplementary 
fig. S9, Supplementary Material online), consistent with 
known native and nonnative alternative initiation codons 
in £ coli (Sacerdot et al. 1996; Sussman et al. 1996). In addition, 
we observed that GUA, GUC, and GUU could serve as weak 
initiation codons (7-14% as efficient AUG). Initiation from 
GUA in £ coli has been previously reported (Haggerty and 
Lovett 1997), but initiation from GUC and GUU has not. 

Mutational Tolerance 

As the effects of nonsynonymous mutations dwarfed that of 
synonymous mutations, we combined the gene fitness data 
of synonymous codons to determine the DFE of missense 
mutations in the TEM-1 protein (fig. 2). This protein fitness 
landscape of TEM-1 broadly matched what is known about 
protein structure in general and TEM-1 in particular. For 
example, proline was the least tolerated substitution (supple- 
mentary fig. S10, Supplementary Material online, which dis- 
plays TEM-1 's amino acid substitution matrix for Amp 
resistance), especially in alpha helices, and key TEM-1 active 
site residues did not tolerate mutation (fig. 2). TEM-Ts signal 
sequence is required for export via the Sec pathway to the 
periplasm. The signal sequence (fig. 2) tolerated most muta- 
tions consistent with the pathway's broad specificity 
(Gierasch 1989). However, the hydrophobic core of the 
signal sequence did not tolerate substitution of charged res- 
idues, consistent with typical export-defective mutants in Sec 
pathway signal sequences (Gierasch 1989). Signal sequence 
residue L21 was a hot spot for beneficial mutations, and L21F 
is found in some extended-spectrum-resistant TEM alleles 
(Sougakoff et al. 1989). 

The comprehensive protein fitness landscape of missense 
mutations enables a rigorous determination of a protein's 
mutational tolerance in its biological context. We determined 
the effective number of amino acids at a position (k*), which 
derives from the fitness entropy that is calculated from the 
distribution of protein fitness values for the 20 amino acids at 
that position. This measure of tolerance is more informative 
than establishing an arbitrary fitness cutoff for deciding 
whether a mutation is tolerated. Our approach is analogous 
to how information-theoretical entropy is used to measure 
variability at a position in a set of aligned sequences (Shenkin 
et al. 1991). However, our measure of tolerance is specific for 



TEM-1 and the effect of single amino acid substitutions. This 
tolerance is a measure of TEM-1 's ability to move a Hamming 
distance of one on the amino acid level and thus does not 
include epistatic effects. A k* value of 1 corresponds to a 
position at which all missense mutations completely inacti- 
vate the protein and a k* value of 20 means that all 19 amino 
acid substitutions provide the same fitness as the wild-type 
amino acid. 

The distribution of k* was strongly biased toward high 
values (fig. 4A). Half of all positions accepted 15.5 or more 
amino acid substitutions. Under the simple assumption of a 
linear correlation, percent solvent-accessible surface area 
accounted for 49% of the k*'s variance (fig. 4B) and predicted 
k* better than distance from the active site (fig. 4C) or a k* 
determined from a sequence alignment of 156 class A 
(3-lactamases (Deng et al. 2012) (fig. 4D). Both a k* based on 
a multiple sequence alignment (fig. 4D) and previous calcula- 
tions of k* for TEM-1 (Deng et al. 2012) (supplementary 
fig. S1 1, Supplementary Material online) greatly underesti- 
mated TEM-1 's mutational tolerance to single amino acid 
substitutions presumably because epistatic constraints will 
further limit what sequence combinations are seen naturally, 
the set of known functional sequences is only a small subset of 
all possible functional sequences, and a high stringency was 
used in selecting functional sequences in the case of the later 
study. The tolerance of amino acid position /' weakly correlated 
with positions /' + 1, i + 3, and i + 4 (correlation coefficient 
0.25-0.28, P < 1.8 x 10~ 5 ) but not / + 2 or /> 4. This corre- 
lation primarily occurred at residues with high k* values. The 
eight positions with k* < 2.5 include the four strictly con- 
served residues involved in the catalytic mechanism (S70, 
K73, S130, and E166) and four other highly conserved residues 
(fig. 4E). In contrast, the 42 most tolerant positions (k* > 19) 
predominantly appeared away from the active site in surface 
loops and at position 2 in alpha helices (fig. 4E). Alpha helices 
(mean k* - 1 3.5 ± 5.4) tolerated substitutions better than beta 
strands (mean k* = 9.89 ±4.8) (P = 0.0005 Student's t-test), 
perhaps a reflection of the buried nature of the beta strands. 

Substitution matrices, such as BLOSUM (Henikoff S and 
Henikoff JG 1992), score the likelihood of substituting one 
amino acid for another based on alignments of conserved 
regions of related proteins. A recent study found that the 
BLOSUM62 matrix best predicted the effect of nonsynon- 
ymous mutations in TEM-1 and explained 16% of the vari- 
ance in the MIC for amoxicillin (Jacquier et al. 2013). We find 
that BLOSUM62 matrix scores predict 16% of the variance in 
protein fitness caused by nonsynonymous mutations. 

Determinants of Mutational Effects on Protein Fitness 
What basic phenomena underlie the DFE? For an enzyme, 
fitness (w) will strongly depend on the total catalytic activity 
in the cell (v t ), which is a product of the enzyme's specific 
catalytic activity (v sp ) and the protein abundance (P), which is 
how much protein is present in the cell in a correctly folded, 
soluble form. 

w=f(v t ) (1) 
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Fic. 4. Tolerance of TEM-1 to missense mutation. Tolerance was measured by the effective number of amino acids at a position (k*), which derives from 
the distribution of protein fitness values for the 20 amino acids at that position, k* ranges in value from 1 (position is completely intolerant to 
substitution) to 20 (position tolerates all possible amino acids with no loss in fitness). (A) The distribution of k* values in TEM-1. (B) Correlation of k* 
correlates with percent solvent-accessible surface. (C) Correlation of k* with distance from the active site. (D) Correlation of k* with a sequence 
alignment of 156 class A (3 lactamases (Deng et al. 2012). (E) Model of TEM-1 (PDB ID 1XPB [Fonze et al. 1995]) indicating the least tolerant positions 
(k* < 2.5, shown in blue), which include the key active site residues S70, K73, S130, and E166, and the most tolerant positions (k* > 19, shown in red). 



Vt = Vs P P (2) 

For many genes, especially essential genes, the functional 
form of equation (1) is likely complex. For example, an 
increase in v t may be detrimental for fitness if it neg- 
atively perturbs metabolic flux in the cell. In addition, essen- 
tial genes are likely to have evolved to be buffered against the 
deleterious effects of mutation (i.e., they possess a v t that is 
well above a level that would compromise fitness). A study of 
27 mutants of the essential enzyme dihydrofolate 
reductase (DHFR) supports this idea (Bershtein et al. 2012). 
The mutations were chosen to have little effect on specific 
catalytic activity but a range of effects on thermostability. The 
authors found that organismal fitness (i.e., growth rate) only 
weakly correlated with protein abundance, and large de- 
creases in protein abundance generally had marginal effects 
on fitness. A follow-up study estimated that DHFR could 
sustain an 80% cut in protein abundance with little effect 
on organismal fitness and that the dependence of fitness 
on protein abundance exhibited Michaelis-Menten-like 
behavior (Bershtein et al. 2013). This study illustrates the 
challenge of gaining insight into the basic phenomena 
underling the DFE without knowledge of the form of 
equation (1). TEM-1 offers a simple case for addressing this 
issue, as TEM-1 fulfills a single cellular role (inactivation of (3- 
lactam antibiotics), and the reaction's substrate and 
product are not part of any native E. coli metabolic or 
signaling pathway. As a result, fitness for TEM-1 is di- 
rectly proportional to the total antibiotic hydrolysis 



activity in the cell (Soskine and Tawfik 2010), as shown in 
equation (3): 

w oc v t = v sp P (3) 

and this facilitates an analysis of the relative effects of muta- 
tion on both v sp and P. 

Experimental evolution studies have shown that v sp and 
P are equally important targets for adaptive evolution 
(Counago et al. 2008; Walkiewicz et al. 2012). We expect 
protein abundance to be a function of the thermodynamic 
stability (AG) as well as protein production rates (i.e., arising 
from mRNA properties, interactions with chaperones) and 
degradation rates (i.e., proteolytic susceptibility). Both com- 
putational and experimental studies show that, on average, 
missense mutations decrease thermodynamic stability 
(Tokuriki et al. 2007). A prevailing hypothesis on the origin 
of deleterious fitness effects of mutation states that thermo- 
dynamic stability is the primary determinant of the DFE 
through its effect on protein abundance (DePristo et al. 
2005; Camps et al. 2007; Tokuriki et al. 2007; Wylie and 
Shakhnovich 2011). Although the hypothesis is intuitive 
and appealing, experimental evidence for a significant corre- 
lation between protein stability and fitness via an effect on 
protein abundance is scant (Bershtein et al. 2012). Mutations 
that reduce function often show decreased protein abun- 
dance (Pakula et al. 1986; Schultz and Richards 1986); how- 
ever, mutations that increase stability can reduce specific 
activity (Shoichet et al. 1995) and reductions in protein 
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Specific catalytic 
activity 

Fic. 5. The determinants of protein fitness. (A) Loss of fitness correlates with loss of thermodynamic stability. Protein fitness is shown as a function of 
change in AC as predicted by Rosetta (Chaudhury et al. 2010) for 4,783 missense mutations in wild-type TEM-1. The median predicted A AC for fitness 
deciles is shown in red triangles. Predicted changes >15 REU are not shown and are not considered in the median calculation. The distribution of A AG 
for select fitness deciles can be found in supplementary figure S12A, Supplementary Material online. (B) Total cellular catalytic activity determines TEM- 
1 protein fitness. The average total cellular catalytic activity <i/ t > and the average protein abundance <P> were experimentally measured for 13 
sublibraries of ~1 5,000 unique TEM-1 alleles partitioned based on relative fitness. The values of <v t > and <P> are relative to that of TEM-1. The slight 
sigmoidal form of <v t > is an expected artifact of the methodology (supplementary fig. S14, Supplementary Material online). The error bars represent 
the standard deviation of six assays from two independent experiments. The lines are guides for the eye. (C) Protein fitness phase space defined by 
equation (3). The protein abundance and specific catalytic activity (relative to TEM-1) of 26 randomly selected members of sublibraries 6 (red solid 
circle) and 7 (blue open square) is shown. The dotted line corresponds to an equal decrease in protein abundance and specific catalytic activity. In 
region 1, a mutation affects specific catalytic activity more than protein abundance. The solid lines are of constant fitness at the average expected fitness 
values of the two sublibraries from which the alleles were randomly selected. 



stability often accompany adaptive mutations (Wang et al. 
2002). In the aforementioned study of 27 mutations that 
destabilized DHFR (most of which were buried in the hydro- 
phobic core of the protein) the authors found that although 
protein abundance correlated with thermostability (^ = 0.41 
at 37 °C), organismal fitness changed very little with protein 
abundance (Bershtein et al. 2012). The study could not 
address how the deleterious effects of mutations partition 
between effects on specific catalytic activity and protein 
abundance because the mutations were selected to be 
those that do not affect catalytic activity. A study of 990 
missense mutations in TEM-1 found that 15-19% of the var- 
iance in amoxicillin resistance could be explained by the com- 
putationally predicted change in protein stability caused by 
the introduction of the mutation in TEM-1 16 (TEM-1 16 is 
TEM-1 with the V84I and A184V mutations) (Jacquier et al. 
2013); however, the study did not address the mutations' 
effect on protein abundance or specific catalytic activity. 
A comprehensive, systematic study of 1) the relationship 
between fitness and thermostability, and 2) the relative con- 
tributions of protein abundance and specific activity to the 
deleterious effects of mutations would more fully address the 
fundamental phenomena underlying the DFE. 

We predicted A AG (AG wi | d . type - AG mutant ) using Rosetta 
(Das and Baker 2008; Chaudhury et al. 2010) for 4,783 mis- 
sense mutations of TEM-1, allowing limited backbone flexi- 
bility (fig. 5A). Variants that were predicted to be more stable 
tended to have higher protein fitness (fig. 5A; supplementary 
fig. S12A, Supplementary Material online). The larger a mu- 
tation's deleterious effect on fitness, the higher the probability 



that the mutation produced a very large predicted energy 
score (supplementary fig. S12A, Supplementary Material on- 
line). Predictions of A AG using PoPMuSiC (Dehouck et al. 
2011), a more empirical approach than Rosetta to predicting 
changes in protein stability, produced similar results and in- 
dicated that 18% of the variance in protein fitness can be 
explained by thermostability (supplementary fig. S12B, 
Supplementary Material online). We compared fitness with 
experimentally measured melting temperatures of 36 TEM-1 
alleles and found a positive correlation with an r 2 of 0.53 
(supplementary fig. S12D, Supplementary Material online). 
This suggests that limitations in the computational prediction 
of A AG result in an underestimation of the degree to which 
thermodynamic stability determines fitness. The lack of a 
positive correlation between melting temperature and fitness 
observed in a previous study (Bershtein et al. 2012) under- 
scores the fact that protein stability effects on fitness will be 
observed only when the fitness function of equation (1) is in a 
regime where changes to v t affect fitness. Our systematic and 
comprehensive approach to examining the relationship be- 
tween protein stability and protein fitness for a single protein, 
combined with our observed correlation between melting 
temperature and fitness provides strong experimental evi- 
dence that effects on protein stability significantly shape 
the DFE. 

Whether mutational reductions in protein abundance as 
opposed to specific activity are the major cause of loss of 
fitness has not previously been experimentally addressed. 
We experimentally addressed this question by analyzing the 
soluble fraction of cell lysates of the 13 sublibraries and 
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randomly selected alleles from our TEM-1 library. We first 
established that w and v t are directly proportional as pre- 
dicted by equation (3) by measuring the mean total hydrolysis 
activity of the cell <v t > for the 13 sublibraries of CCM-2 
(fig. 5B). We measured protein abundance by quantitative 
western blot of the soluble fraction of cell lysates (supplemen- 
tary fig. S13, Supplementary Material online). We assumed all 
soluble protein was folded and active. The mean protein 
abundance <P> of the sublibraries did not decrease nearly 
as rapidly with decreasing protein fitness as <i/ t > did (fig. 5B). 
In addition, an increase in the mean amount of aggregated 
protein did not accompany a loss of fitness (supplementary 
fig. S15, Supplementary Material online). These findings sug- 
gest that mutational effects on v sp rather than on P may play 
the larger role in the deleterious effects of mutation. As this 
interpretation hinges on the distribution of values of P in the 
sublibraries, we measured P for 27 randomly selected alleles 
with a protein fitness of about 0.025-0.05 (i.e., the alleles of 
supplementary fig. S4, Supplementary Material online). We 
chose this fitness range so that the mutational effects were 
substantial, but not inactivating. This ensured that our con- 
clusions would not depend on small changes in w and P. We 
excluded the I13E allele from analysis since this mutation 
in the signal sequence caused a defect in normal export/pro- 
teolytic processing (supplementary fig. S13, Supplementary 
Material online). The remaining 26 alleles exhibited a decrease 
in both protein abundance and predicted thermodynamic 
stability relative to TEM-1 with the exception of the R244E 
allele, which showed an increase in both (supplementary 
fig. S12C, Supplementary Material online). From w and P, 
we calculated v sp using equation (3) and examined the pro- 
tein fitness phase space by plotting P versus v sp (fig. 5C). We 
find that the deleterious effects of mutations, at least for 
mutations with large deleterious effects, arise more from a 
decrease in specific activity than from a decrease in protein 
abundance. Despite the large negative effects on specific ac- 
tivity, the mutated residues of the 26 alleles were not clus- 
tered around the active site but were scattered throughout 
the protein (supplementary fig. S4C, Supplementary Material 
online). Thus, the dominant effect of mutation on specific 
activity does not arise because the 26 mutations were biased 
to be proximal to the active site. We postulate that muta- 
tional effects on specific activity may be as important to the 
DFE at high fitness as at low fitness, but this postulate requires 
experimental investigation. 

We do not interpret the diminished role of mutational 
effects on protein abundance as reducing the role of thermo- 
dynamic stability in fitness. Protein stability, in addition to its 
effect on protein abundance, may exert its effect on fitness 
through a decrease in a protein's specific activity. Perhaps, this 
manifests by perturbing the conformational ensemble away 
from more active states or by increasing the number of states 
(i.e., altering protein dynamics). Protein abundance's relative 
resilience to decreases in thermodynamic stability is striking 
but fits the growing appreciation that the cellular environ- 
ment is not a passive solution at equilibrium (Bershtein et al. 
2013). Rather, the cellular environment acts as a buffer for 
deleterious mutational effects on protein abundance through 



the effect of chaperones, proteins that facilitate the proper 
folding of other proteins (Rutherford 2003). Chaperone over- 
expression can compensate for the deleterious mutational 
effects on protein abundance (Tokuriki and Tawfik 2009) 
and fitness (Bershtein et al. 2013). This theory offers an 
explanation for TEM-1 's stability threshold that buffers the 
effect of mutations on fitness (Bershtein et al. 2006). We 
suspect that the relative contribution of protein abundance 
to fitness may increase with the number of mutations as the 
protein's stability margin is exhausted by the cumulative 
effect of mutations, an effect that is characterized by negative 
epistasis (Bershtein et al. 2006). As such, negative epistasis 
may arise in part as a consequence of the beneficial properties 
of the cellular environment in addition to a protein's intrinsic 
stability margin. 

Conclusions 

The application of synthetic biology to the study of funda- 
mental biological questions, as we have done in this study of 
gene and protein fitness landscapes, offers a well-defined, 
systematic approach for testing and generating hypotheses. 
Our comprehensive determination of the fitness effects of 
mutation of TEM-1 provides the first detailed maps of fitness 
landscapes corresponding to a gene and its nearest neighbors 
at the basepair, codon, and amino acid level. To the extent 
that TEM-1 is a representative gene, our study provides several 
important insights. Evolution must traverse fitness landscapes 
under the constraints of the genetic code — constraints that 
minimize the effect of mutation and enrich for beneficial 
mutations. The small fitness effects of synonymous mutations 
have complex determinants including regional proclivities for 
synonymous fitness effects in the gene. At the beginning of 
the gene, fitness effects of synonymous mutations strongly 
correlate with mRNA stability. Missense mutational effects on 
thermodynamic stability shape the DFE, but their deleterious 
effect on specific protein activity tends to exceed that on 
protein abundance, at least for mutations with large delete- 
rious effects. We hypothesize that TEM-1's high mutational 
tolerance may in part derive from the cell's buffering capacity 
to mediate the deleterious effects of lost stability on protein 
abundance, a phenomena that might give rise to negative 
epistasis. Further inquiry into the fundamental determinants 
of the landscape's topology is necessary to address this 
hypothesis and substantiate these findings. 

Materials and Methods 

Fitness Determination 

Escherichia coli SNO301 (ampDl, ampAl, ampC8, pyrB, reck, 
and rpsL) cells harboring the comprehensive codon mutagen- 
esis library CCM2 (Firnberg and Ostermeier 2012) were plated 
on LB-agar plates supplemented with 13 different Amp con- 
centrations (2-fold increments ranging from 0.25ug/ml to 
1,024 ug/ml) to partition the library into 13 partially over- 
lapping sublibraries based on relative Amp resistance using 
a synthetic gene circuit that functions as a tunable band-pass 
genetic selection for Amp resistance (Sohka et al. 2009) (sup- 
plementary figs. SI and S2, Supplementary Material online). 
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Barcoded PCR amplicons were prepared from each plate and 
subjected to 454 deep-sequencing. The 1,325,979 sequencing 
reads were filtered for quality and for reads that only con- 
tained one codon substitution. We tabulated the number of 
sequencing counts for each allele at each Amp concentration 
and determined the fitness w relative to TEA/1-7 from the 
statistics. As the distribution of growth as a function of 
Amp is roughly symmetric when plotted as the log 2 (Amp 
concentration) (Sohka et al. 2009), we determined the unnor- 
malized fitness/of allele / as 



E c ',p |o g 2 ( a P ) 
p=1 

13 

E c i, P 



(4) 



in which c liP is the number of counts of allele / on sublibrary 
plate p in the deep sequencing data and a p is the concentra- 
tion of Amp on sublibrary plate p in (.ig/ml. We normalized all 
fitnesses by the fitness of wild type as follows: 



2/wT 



(5) 



This result is a normalized fitness w-, that is 1.0 for wild-type 
TEM-1, greater than 1.0 for beneficial mutations, and between 
0 and 1.0 for deleterious mutations. We determined the fit- 
ness of wild-type TEM-1 (fwr) using equation (4) using the 
counts of all alleles with a synonymous substitution in TEM-1, 
because the fitness of these varied very little. As a check, we 
compared this value with the fitness determined by equation 
(4) using the counts of all sequencing reads that lacked a 
mutation. The two values differed by only 2.5%. We deter- 
mined an upper limit on the error in our fitness measure- 
ments from the DFEs of synonymous mutations as a function 
of the number of times an allele was observed (supplemen- 
tary fig. S3, Supplementary Material online). Fitness values and 
error estimates are tabulated in supplementary data S1 and 
S2 (Supplementary Material online). 

Prediction of mRNA Stability at the Transcript Start 
The RNAfold utility of the Vienna RNA Package (version 2.1.2) 
was used to predict the minimum free energy of RNA se- 
quences (Hofacker 2009). For each allele in supplementary 
figure S7 (Supplementary Material online), the Gibbs free 
energy was calculated as the average free energy of every 
39 nt window centered on nucleotides from —5 to + 10 of 
the gene start as described (Bentele et al. 2013). 

Mutational Tolerance 

The observed effective number of amino acids /c* at a position 
was determined from the protein fitness values of the n mis- 
sense mutations with fitness data at that position using equa- 
tions (6) and (7). 



K = 2 00 

We obtained the effective number of amino acids k* by nor- 
malizing k* to be based on 20 amino acids by equation (8). 



20k* 



(8) 



w, log 2 



(6) 



A table of k* a and k* is provided as supplementary data S4 
(Supplementary Material online). 

Prediction of Protein Thermodynamic Stability 
PyRosetta v3.4.0 r55307 (Chaudhury et al. 2010) was used to 
compute the difference in score (in Rosetta energy units 
[REU]) between the mature structures (lacking the signal se- 
quence) of each amino acid mutant and wild-type TEM-1 
(Protein Data Bank identifier 1XPB; Fonze et al. 1995). 
PopMusic predictions of A AG (supplementary fig. S12B, 
Supplementary Material online) were determine online at 
http://babylone.ulb.ac.be/popmusic (last accessed February 
17, 2014) (Dehouck et al. 2011) using 1XPB. 

Protein Abundance and Total Catalytic Activity 
Relative protein abundance was quantified by using Quantity 
One 1-D analysis software (Bio-Rad) of Western blots of the 
soluble fraction of cell lysates in comparison with a standard 
curve. Representative westerns are shown in supplementary 
figure S13, Supplementary Material online. Catalytic activity 
of the sublibraries and clones was determined by measuring 
the rate of hydrolysis (at 486 nm) of 50 uM nitrocefin in 
10 mM phosphate buffer pH 7.4 at 37 °C The initial rate 
was normalized by the total amount of protein added for 
each sample. 

Supplementary Material 

Supplementary materials and methods, figures SI— SI 5, and 
data S1-S4 are available at Molecular Biology and Evolution 
online (http://www.mbe.oxfordjournals.org/). 
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